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HIGH-ORDER  DISCONTINUOUS  GALERKIN  METHODS: 
SIMULATION  OF  COIL  FLOWS 


I.  PIVKIN,  R.M.  KIRBY  AND  G.E.  KARNIADAKIS 
Division  of  Applied  Mathematics 
Brown  University 


Abstract. 

We  have  developed  spectral  methods  in  the  discontinous  Galerkin  frame- 
work appropriate  for  simulations  of  high-speed  flows  in  complex-geometry 
domains.  In  this  paper  we  present  details  of  the  stability  of  the  method 
and  demonstrate  the  importance  of  over-integration  for  strongly  nonlinear 
problems.  We  then  present  results  from  the  application  of  the  method  to 
stability  studies  of  supersonic  and  subsonic  flows  in  a Chemical  Oxygen 
Iodine  Laser  (COIL)  device. 


1.  Introduction 

Finite  volume  methods  have  been  very  successful  in  simulating  steady 
high-speed  flows  but  they  arc  rather  inefficient  for  unsteady  flow  simula- 
tions, and  especially  for  direct  or  large  eddy  simulations  of  turbulent  and 
transitional  compressible  flows.  Discontinuous  Galerkin  Methods  (DGM) 
[1],  when  combined  with  high-order  discretizations  as  in  [2],  offer  some  of 
the  advantages  of  finite  volume  methods  and  lead  to  numerical  solutions 
with  significantly  reduced  numerical  dispersion  and  dissipation.  In  sum- 
mary, high-order  DGM  are: 

— High-order  finite  volumes  applied  to  structured  or  unstructured  meshes. 

— Flux-based,  and  thus  maintain  conservativity  which  is  important  for 
correct  shock  location  and  long-time  integration. 

— L'2-stable,  and  thus  they  do  not  require  explicit  flux  limiters. 

— Robust  as  they  employ  Rieinann  solvers. 

The  particular  DGM  implementation  we  have  developed  employs  an 
orthogonal  polynomial  basis  of  different  order  in  each  element.  The  discon- 
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tinuous  basis  is  orthogonal,  hierarchical,  and  maintains  a tensor-product 
property  even  for  non-separable  domains  [3,  4]. 

Unlike  pseudo-spectral  (collocation)  methods  used  often  for  incompress- 
ible flow  simulations,  in  our  spectral  DGM  approach  the  number  of  quadra- 
ture points  and  the  number  of  degrees  of  freedom  (i.e.,  modes)  are  decoupled 
in  each  element.  Specifically,  it  is  a super-collocation  method  combined  with 
Galerkin  projections  in  an  element-wise  fashion.  This  is  very  important  for 
the  treatment  of  nonlinear  terms  as  this  approach  offers  the  possibility  of 
dealiasing  on  arbitrarily  non-uniform  grids.  We  discuss  this  issue  in  detail 
in  the  next  section,  where  we  show  that  simple  quadrature-modes  rules  can 
be  derived  by  considering  the  long-term  (asymptotic)  stability  of  DGM. 

In  the  last  section,  we  apply  DGM  to  three-dimensional  supersonic  and 
subsonic  flows  in  a Chemical  Oxygen  Iodine  Laser  (COIL)  configuration. 
The  flow  field  of  a COIL  typically  contains  multi-phase  flow  (oxygen,  io- 
dine and  helium)  as  well  as  chemical  reactions  [5] , but  here  we  consider  the 
stability  of  cold  (helium)  flow.  Due  to  the  importance  of  COILs  in  military 
and  industrial  uses,  there  has  been  a lot  of  research  focused  on  them  during 
the  last  decade  [6,  7,  8].  Recent  two-dimensional  simulations  have  shown 
that  the  supersonic  COIL  flow  is  unsteady  although  corresponding  pre- 
liminary simulations  have  shown  substantially  reduced  temporal  variations 
[9].  In  the  current  work,  we  consider  both  supersonic  as  well  as  subsonic 
conditions  on  different  grids  and  with  p-refinement  to  address  more  sys- 
tematically the  question  of  unsteadiness.  Our  results  show  that  the  strong 
compressibility  and  the  symmetric  crossflow  stabilize  the  supersonic  nozzle 
flow,  in  contrast  with  the  subsonic  cases. 

2.  Discontinuous  Galerkin  Method  - DGM 

We  consider  the  non-dimensional  compressible  Navier-Stokes  equations, 
which  we  write  in  compact  form  as 

17t  -1-  V • F = Re^V  ■ F"  in  D (1) 

where  F and  F"  correspond  to  inviscid  and  viscous  flux  contributions,  re- 
spectively, and  Reoo  is  the  reference  Reynolds  number.  Here  the  vector 
U = [p, pui, pu2, pus, pe]1  with  u = (ui,U2,U3)  the  local  fluid  velocity,  p 
the  fluid  density,  and  e the  total  energy.  Splitting  the  Navier-Stokes  opera- 
tor in  this  form  allows  for  a separate  treatment  of  the  inviscid  and  viscous 
contributions,  which,  in  general,  exhibit  different  mathematical  properties. 

To  give  an  overview  of  the  formulation  we  first  apply  DGM  to  the  linear 
two-dimensional  equation  for  advection  of  a conserved  quantity  u 


+ V-F(u)  = 0, 


(2) 
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where  F(u)  = ( f(u),g(u),h(u ))  is  the  flux  vector  which  defines  the  trans- 
port of  u(x,  t).  We  start  with  the  variational  statement  of  the  standard 
Galerkin  formulation  of  (2)  by  multiplying  by  a test  function  v and  inte- 
grating by  parts 

[ ^-vdx  + [ vh-F(u)ds—  [ Vn-F(it)  dx  = 0.  (3) 

Jo.  ot  Jen  Jn 

The  solution  u € X (approximation  space)  satisfies  this  equation  for  all  v 6 
V (test  space),  where  X may  contain  discontinuous  functions.  The  discrete 
space  X6  contains  polynomials  within  each  “element,”  but  zero  outside 
the  element.  Here  the  “element”  is,  for  example,  an  individual  triangular 
region  Tj  in  the  computational  mesh  applied  to  the  problem.  Thus,  the 
computational  domain  0 = (J?:  Tj,  and  Tj,Tj  overlap  only  on  edges. 


Figure  1.  Interface  conditions  between  two  adjacent  triangles. 

Each  element  ( E ) is  treated  separately,  giving  a variational  statement 
(after  integrating  by  parts  once  more): 

J^( u,v)e  + J ^ v{f{ui,ue)  - F{ui))-nds  + (V  • F{u),v)E  = 0,  (4) 

where  F(ui)  is  the  flux  of  the  interior  values.  Computations  on  each  el- 
ement are  performed  separately,  and  the  connection  between  elements  is 
a result  of  the  way  boundary  conditions  are  applied.  Here,  boundary  con- 
ditions are  enforced  via  the  numerical  surface  flux  f(ui,ue ) that  appears 
in  equation  (4).  Because  this  value  is  computed  at  the  boundary  between 
adjacent  elements,  it  may  be  computed  from  the  value  of  u given  at  either 
element.  These  two  possible  values  are  denoted  here  as  Ui  in  the  interior  of 
the  element  under  consideration  and  and  ue  in  the  exterior  (see  figure  1). 
Upwinding  considerations  dictate  how  this  flux  is  computed. 

In  the  more  complicated  case  of  a hyperbolic  system  of  equations,  an 
approximate  Riemann  solver  should  be  used  to  compute  a value  of  /,  g,  h 
(in  three-dimensions)  based  on  w*  and  ue.  Specifically,  we  compute  the  flux 
f(ui,ue ) using  upwinding,  i.e. 

f(u)  = RK^Lui  + Rh~  Lue 
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where  A (the  Jacobian  matrix  of  F)  is  written  in  terms  of  the  left  and  right 
eigenvectors,  i.e.  A = RAL  with  A containing  the  corresponding  eigenvalues 
in  the  diagonal;  also,  A±  = (A±|A|)/2.  Alternatively,  we  can  use  a standard 
Lax-Friedrichs  flux 

f(u)  = ^(/K)  + /(lit))  - ^R\A\L(ue  - Ui). 

In  order  to  solve  the  compressible  Navier-Stokes  equations  we  also  need 
to  treat  the  second-order  terms.  This  is  done  similarly  by  introducing  aux- 
iliary fluxes  and  constructing  a system  of  first-order  equations  that  we  treat 
similarly  as  the  advection  equation  above.  The  difference  is  in  the  choice 
of  the  numerical  fluxes  at  the  interface,  which  we  have  taken  here  to  be 
the  arithmetic  mean  of  the  fluxes  from  adjoint  elements.  However,  other 
choices  are  possible  that  can  affect  the  computational  complexity  and  the 
accuracy  [1]. 

With  regard  to  the  trial  basis  we  employ  Jacobi  polynomials  with  vari- 
able weights  as  reported  in  [2].  The  degree  of  polynomial  can  be  varied 
from  one  element  to  the  other  but  interpolation  along  edges  or  elemental 
faces  should  involve  the  higher  order  polynomial  in  order  to  maintain  con- 
servativity.  The  DGM  method  is  I^-stable  but  in  practice  instabilities  may 
arise  in  under-resolved  or  marginally  resolved  simulations.  Such  instabili- 
ties were  first  reported  in  [2]  where  a slight  over-integration  seem  to  provide 
long-time  integration  stability.  In  the  next  section,  we  revisit  this  issue  and 
provide  more  details  on  over-integration  for  both  quadratic  as  well  as  cubic 
nonlinearities. 

2.1.  STABILITY  AND  OVER-INTEGRATION 

To  understand  the  ramifications  of  under-integration  of  nonlinear  terms, 
we  performed  the  following  test: 

1.  Initialize  a single  element  spanning  [—1,1]  and  containing  16  modes. 

2.  Initialize  all  the  modal  coefficients  to  one. 

3.  Evaluate  the  modal  representation  on  a set  of  quadrature  points  q. 

4.  Pointwise  square  the  values  at  the  quadrature  points. 

5.  Pre-multiply  the  set  of  points  (as  a vector)  by  the  collocation  derivative 
matrix  of  the  appropriate  size  (rank  q x q). 

6.  Project  back  to  modal  coefficients  by  discrete  inner  products  using 
Gaussian  integration. 

The  procedure  above  mimicks  the  “physical  space”  or  pseudo-spectral 
evaluation  of  the  term  ^ commonly  used  in  spectral  methods  for  evaluat- 
ing nonlinear  terms.  This  test  was  chosen  because  even  in  its  simplicity  it 
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models  the  order  of  nonlinearity  that  occurs  in  the  solution  of  the  incom- 
pressible Navier-Stokes  equations.  All  modes  are  set  to  one  to  mimick  a case 
in  which  an  element  has  under-resolved  or  marginally  resolved  the  solution 
within  the  element.  In  the  test  above,  the  only  unspecified  parameter  is  the 
number  of  quadrature  points  q to  be  used.  In  using  Gauss-Lobattto  points, 
the  value  of  q is  taken  to  be  one  more  than  the  number  of  modes  M (in 
this  case  then  M = 16  and  q = 17)  [2],  but  this  value  is  appropriate  for 
the  inner  products  corresponding  to  linear  terms.  For  a quadratic  or  cu- 
bic nonlinearities  more  quadrature  points  are  required.  The  ramifications 
of  under-integration  of  this  form  are  shown  in  figure  2.  The  figure  on  the 
left  was  obtained  for  quadratic  nonlinearity  and  the  figure  on  the 

right  was  obtained  for  a cubic  nonlinearity  ( The  difference  in  the 
modal  coefficients  at  the  conclusion  of  the  algorithm  above  for  different  val- 
ues of  q is  provided.  We  observe  that  for  the  quadratic  nonlinearity,  once 
| M quadrature  points  are  used,  the  difference  in  the  modal  values  do  not 
change.  Similarly  for  the  cubic  nonlinearity,  once  2 M quadrature  points  are 
used,  the  difference  in  the  modal  values  do  not  change. 
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Figure  2.  Comparison  of  the  difference  in  modal  coefficients  when  different  numbers  of 
quadrature  points  are  used.  Quadratic  nonlinearity  on  the  left  and  cubic  nonlinearity  on 
the  right. 

If  figure  3 we  plot  the  difference  between  using  the  |M  rule  and  2 M 
rule  (left  and  right  figures,  respectively)  versus  M + 2 rule  for  chosing  the 
number  of  quadrature  points  q.  Observe  that  for  small  number  of  modes 
the  two  regions  overlap  or  may  be  sufficiently  close  that  using  M + 2 will 
not  lead  to  aliasing  instabilities.  This  may  be  an  explanation  of  the  results 
shown  in  [2]  in  which  over-integration  by  one  point  was  sufficient  to  stabilize 
the  flow  simulation. 

To  further  test  the  integration  of  the  nonlinear  terms,  we  chose  to  solve 
viscous  Burgers  equation  with  v = 10-5.  Five  equally  spaced  elements 
spanning  [—1,1]  were  used,  each  one  having  16  modes.  We  have  found 
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Figure  3.  Difference  between  using  the  | M rule  and  2 M rule  (left  and  right  figures, 
respectively)  versus  (M  + 2)  rule  for  chosing  the  number  of  quadrature  points  q. 


that  when  using  17,  19  and  21  quadrature  points,  the  solution  is  unstable 
(measured  by  the  L 2 norm).  Once  the  number  of  quadrature  points  reaches 
24  (§M,  where  M is  the  number  of  modes),  the  L2  norm  of  the  solution 
does  not  change. 

3.  Simulation  of  COIL  Flows 

The  chemical  oxygen  iodine  laser  (COIL)  is  a very  powerful  laser,  capable  of 
producing  megawatts  of  continuous  power  at  short  wavelengths  (1315  nm). 
There  are  two  distinct  types  of  COIL  configurations,  depending  on  the 
characteristic  flow  velocity  of  the  constituent  gases.  In  the  first  one,  the 
supersonic  lasers,  the  velocities  are  of  the  order  of  400  m/s  or  higher  while 
in  the  second  type,  the  subsonic  lasers,  the  characteristic  velocities  are 
an  order  of  magnitude  lower.  There  has  been  some  uncertainty  regarding 
the  stability  of  the  COIL  flows,  especially  in  the  supersonic  regime.  Two- 
dimensional  computations  show  that  the  flow  is  unstable  to  small  perturba- 
tions and  becomes  oscillatory  with  frequency  of  about  40  Khz  [9].  On  the 
other  hand,  preliminary  three-dimensional  simulations  reported  recently  in 
[9]  show  only  extremely  small  time  variations  in  amplitude.  Unfortunately, 
there  are  not  enough  experimental  measurements  or  flow  visualizations  to 
document  the  stability  of  the  COIL  flows,  especially  in  the  supersonic  noz- 
zle (the  so-called  RADICL  nozzle).  In  this  work,  we  study  the  stability  of 
COIL  flows  via  direct  numerical  simulations  based  on  the  DGM  spectral/fip 
element  method  and  the  code  NEKTAR  described  in  the  previous  section; 
see  also  [1,  4]. 

An  overview  of  the  supersonic  flow  in  the  central  part  of  the  nozzle  is 
shown  in  figure  4.  There  is  a strong  intercation  between  the  incoming  al- 
most uniform  flow  and  the  crossflow  emanating  from  the  six  nozzles.  The 


COIL  FLOWS 


151 


SUPERSONIC  (CASE  A) 


Inflow 

Pipes 

Outflow 

Density 

0.013066498 

Density 

0.039323356 

Density 

0.001355714 

U-velo 

153.3355557 

U-velo 

442.6130539 

U-velo 

1208.084843 

V-velo 

0 

V-velo 

0 

V-velo 

0 

W-velo 

0 

W-velo 

0 

W-velo 

0 

Energy 

12209.51111 

Energy 

53746.82174 

Energy 

1158.745548 

Mach 

0.152511337 

Mach 

0.38 

Mach 

3.264908 

SUBSONIC  (NO  PIPES)  (CASE  B) 


Inflow 

Pipes 

Outflow 

Density 

0.013066498 

Density 

0 

Density 

0.0130896 

U-velo 

153.3355557 

U-velo 

0 

U-velo 

121.480094 

V-velo 

0 

V-velo 

0 

V-velo 

0 

W-velo 

0 

W-velo 

0 

W-velo 

0 

Energy 

29697.27213 

Energy 

0 

Energy 

29727.77876 

Mach 

0.097424 

Mach 

0 

Mach 

0.077138 

SUBSONIC  (PIPES)  (CASE  C) 


Inflow 

Pipes 

Outflow 

Density 

0.013066498 

Density 

0.039323356 

Density 

0.013445426 

U-velo 

153.3355557 

U-velo 

0 

U-velo 

121.480094 

V-velo 

0 

W-velo 

50.95079845 

V-velo 

0 

W-velo 

0 

W-velo 

0 

W-velo 

0 

Energy 

29697.27213 

Energy 

119387.5174 

Energy 

30553.32435 

Mach 

0.097424 

Mach 

0.2427434 

Mach 

0.077138 

INCOMPRESSIBLE  (CASE  D) 

Inflow 

Pipes 

Outflow 

U-velo 

2.003558732 

U-velo 

0 

U-velo 

0 

V-velo 

0 

V-velo 

17.40503091 

V-velo 

0 

W-velo 

0 

W-velo 

0 

W-velo 

0 

geometry  and  dimensions  of  the  entire  device  are  shown  in  figure  5.  This 
three-dimensional  “slice”  is  repeated  twenty  times  (in  the  direction  per- 
pendicular to  the  page)  in  order  to  make  the  entire  device,  so  there  are 
20  large  diameter  nozzles  and  40  small  diameter  nozzles  in  the  COIL  de- 
vice. We  simulate  only  one  slice,  as  shown  in  the  figure,  and  apply  periodic 
boundary  conditions  along  the  third  direction.  A typical  mesh  employed 
in  the  supersonic  case  is  shown  in  figure  6.  It  consists  of  10,224  hexahedra 
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elements  with  8 layers  in  the  z-  (periodic)  direction  of  variable  order  as 
shown  in  the  figure.  The  subsonic  cases  were  simulated  with  similar  meshes 
but  with  7,312  hexahedra  of  constant  polynomial  order  p = 3.  The  incom- 
pressible case  was  simulated  on  a mesh  consisted  of  3, 066  tetrahedra  with 
polynomial  order  p = 6. 

The  specific  flow  conditions  we  consider  are  described  in  table  3 (all 
units  in  metric  system).  We  have  simulated  in  detail  the  following  four 
cases:  (A)  Supersonic  flow  corresponding  to  experimental  case  described  as 
“9257cfll”  in  experiments  with  helium;  (B)  Subsonic  flow  without  trans- 
verse (secondary)  flow;  (C)  Subsonic  flow;  and  (D)  Incompressible  flow.  In 
all  cases,  the  geometry  remains  identical  and  so  is  the  primary  incoming 
flow.  Only  the  energy  input  and  mass  flowrate  from  the  side  pipes  varies 
in  cases  (2)  and  (3),  corresponding  to  zero  and  approximately  one-eighth 
of  the  primary  flow,  respectively.  In  cases  (A)  and  (D)  the  secondary  flow 
is  approximately  one-fourth  of  the  primary  flow. 


Figure  4 • Streamlines  in  the  central  portion  of  the  RADICL  supersonic  nozzle. 

Based  on  systematic  direct  numerical  simulations,  we  have  concluded 
that  unsteadiness  is  suppressed  in  the  supersonic  conditions  similar  to  the 
ones  employed  in  the  experiments.  A typical  plot  of  Mach  contours  is  shown 
in  figure  7.  Pressure  distribution  along  the  wall  of  the  COIL  as  well  as  along 
the  centerline  of  the  COIL  are  shown  in  figure  8.  Clearly,  there  is  very  good 
agreement  of  the  high-order  DGM  results  with  available  experimental  data 
for  the  wall  pressure  in  contrast  with  corresponding  finite  volume  simula- 
tions performed  here  on  the  same  mesh  of  figure  6 but  with  p = 1.  These 
results  were  obtained  from  converged  (in-time)  simulations  of  the  super- 
sonic nozzle.  To  further  examine  the  stability  of  the  supersonic  flow  we  in- 
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I All  dimensions  in  centimeters  I / Pipes  \ 

1 / (view  from  the  top)  \ 

/ \ 


Figure  5.  Geometry  and  dimensions  of  the  COIL  device. 


Figure  6.  Mesh  showing  the  macro-elements  and  corresponding  polynomial  order  in  the 
central  part  of  the  COIL  (supersonic  case). 


troduced  an  abrupt  and  rather  large  perturbation  into  the  flow  in  the  form 
of  forces  in  all  three  directions  corresponding  to  —8,322  m/s2,  8,322  m/s2 
and  8,322  m/s2  along  the  x-(streamwise),  y-(crossflow)  and  z-(periodic) 
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Figure  7.  Mach  contours  in  the  supersonic  EADICL  nozzle. 


Figure  8.  Pressure  distribution  along  the  wall  and  along  the  centerline  for  the  supersonic 
case. 


directions,  respectively.  This  forcing  was  imposed  for  2.39  x 10~4  ms  and 
subsequently  it  was  removed  while  the  simulation  continued.  A typical  re- 
sult of  the  response  in  the  y-momentum  is  shown  in  figure  9.  We  see  that 
indeed  the  flow  returns  to  steady  state  within  a short  time  interval. 

Unlike  the  supersonic  case,  both  subsonic  cases  we  simulated  converged 
to  time-dependent  flows  while  the  incompressible  flow  transitioned  to  a tur- 
bulent state.  A typical  result  is  shown  in  figure  10  that  plots  Mach  contours 
of  the  subsonic  case  with  crossflow;  a large  scale  unsteadiness  is  present. 
This  is  shown  more  clearly  in  figure  11  where  we  plot  the  time-histories  of 
the  y-momentum  of  both  subsonic  cases.  The  case  without  crossflow  shows 
a stationary  time-periodic  flow  whereas  the  subsonic  flow  with  crossflow 
exhibits  an  additional  modulation  associated  with  the  large  scale  unsteadi- 
ness. Therefore,  it  seems  that  the  effect  of  crossflow  is  to  suppress  unsteadi- 
ness. 
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Figure  9.  Time  history  of  y-momentum  for  the  supersonic  nozzle. 
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Figure  10.  Mach  contours  for  the  subsonic  case  with  crossflow. 
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Figure  11.  Time-history  of  y-momentum  for  the  two  subsonic  cases,  with  crossflow 
(left),  and  without  crossflow  (right). 
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Finally,  we  conclude  by  commenting  on  the  computational  cost  of  the 
simulations.  All  simulations  were  run  using  an  MPI-based  parallel  version 
of  the  method  presented  here  with  the  partitioning  based  on  a multi-level 
graph  approach  provided  by  the  METIS  software  [10].  Specifically,  the  su- 
personic simulation  was  run  on  the  IBM  Power  3 with  0.6  seconds  per  time 
step  on  80  processors,  the  subsonic  simulations  were  run  on  the  IBM  Power 
3 with  0.2  seconds  per  time  step  on  80  processors,  and  the  incompressible 
simulation  was  run  on  the  IBM  SP  604e  (silver  nodes)  with  10.7  seconds 
per  time  step  on  64  processors. 
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